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ABSTRACT 

We  introduce  a new  variance  reduction  technique,  called 
internal  control  variables,  to  be  used  in  the  context  of  regeneration 
simulations.  The  idea  is  to  identify  a sequence  of  control  random 
variables,  each  one  defined  within  a regenerative  cycle,  whose  mean 
can  be  calculated  analytically.  These  controls  should  be  highly 
correlated  with  the  usual  quantities  observed  in  a regenerative 
simulation.  This  correlation  reduces  the  variance  of  the  estimate 


for  the  parameter  of  interest.  Numerical  examples  are  included  for 
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REGENERATIVE  SIMULATION  WITH  INTERNAL  CONTROLS 

by 

Donald  L.  Iglehart  and  Peter  A.  W.  Lewis 

1 . Introduction  and  Suimnary 

Simulators  are  frequently  faced  with  the  task  of  esti- 
mating a parameter  associated  with  the  limiting  distribution  of 
a stable  stochastic  process  which  is  being  simulated.  A 

methodology  based  on  regenerative  processes  for  obtaining  point 
estimates  and  confidence  intervals  for  such  parameters  from  a 
single  realization  of  the  process  has  recently  been  developed 
in  Crane  and  Iglehart  (1974a, b) , (1975a, b)  and  Iglehart  (1975), 

(1976a, b) , and  (1977) . In  this  paper  we  shall  introduce  a new 
technique,  internal  control  variables,  which  can  be  used  in  con- 
junction with  the  regenerative  method  for  obtaining  additional 
variance  reduction  for  the  estimates . 

Suppose  regenerative  process  {X^:t  > 0},  which  is  stable 
in  the  sense  that  = X as  t (here  =»  denotes  weak 

convergence) , is  being  simulated  by  generating  a single  realization 
of  the  process.  For  convenience  think  of  this  process  as  being 
either  a positive  recurrent  Markov  chain  or  the  waiting  time 
process  for  a single  server  queue  with  traffic  intensity  less 
than  one.  Then  simulators  are  freauently  interested  in  estimating 
r = E{f(X)},  for  a given  function  f.  The  principal  goal  of  the 
regenerative  method  is  to  produce  a confidence  interval  estimate 
of  r.  The  regenerative  method  begins  by  observing  the  pairs 
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is  the 


of  random  variables  { :1  £ k £ n}  , where  Tj^ 

length  of  the  kth  regenerative  cycle  and  Yj^  is  the  area  under 
the  function  f (X^)  in  the  kth  cycle.  Two  basic  facts  are  crucial 
for  the  regenerative  method  of  simulating  stable  processes.  First, 
the  pairs  {(Yj^,Tj^):1  <.  k ^ n}  are  independent  and  identically 
distributed  (i.i.d.),  and  second  r = E{f(X)}  = E{Yj^}/E{tj^}  . 

This  last  fact  suggests  that  a natural,  strongly  consistent  point 
estimate  for  r is  r(n)  = y(n)/T{n),  where  Y(n)  = n ^ ^k=l  ^k 
and  T (n)  = n ^ ^k=l^k‘  form  a confidence  interval  for  r 

we  can  use  the  central  limit  theorem  (c.l.t.) 


(1.1) 


/n  (r(n)  - r)  / (a/E{T  ^ N(0,1)  , 


where  N(0,1)  denotes  a mean  zero,  variance  one  normal  random  variable 
2 2 

and  a = E{  [Y^^  - }.  Of  course,  the  constant  a/E{Tj^}  will 

have  to  be  estimated  in  most  simulations. 

The  idea  behind  internal  control  variables  is  to  introduce 
a third  sequence  of  i.i.d.  random  variables  S.  ^ i.  which 

has  the  property  that  Cj^  is  defined  in  terms  of  the  kth  cycle 
and  known  (or  can  be  calculated  analytically)  . Then 

another  strongly  consistent  estimate  for  r is 


I"  1 

(1.2)  r^^(n)  = ^ 

The  subscript  CT  is  meant  to  denote 
ratio  estimator;  similar  estimates 


[Yj^  + - E{Cj^})] 

T (n) 

"controlling  the  top 
will  be  defined  for  " 


" in  the 
bottom  control . 
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with  a replaced 


} 


i 

I 


r 


A c.l.t.  analogous  to  (1.1)  also  holds  for  r^^ 
by  a',  say.  Since  we  are  still  free  to  select  3,  we  choose  to 
do  so  in  such  a way  as  to  minimize  o'.  Having  done  that  we  find 
that 

(o')"  = o^[l  - p^(C,  ,Y,  - rxj]  . 


To  obtain  significant  variance  reduction  a control,  must  be 

found  which  is  highly  correlated  with  a task  that  is 

not  always  easy. 

Section  2 of  the  paper  develops  the  method  of  internal  control 
variables  in  detail  and  contains  specific  examples  associated  with 
the  single  server  queue  and  Markov  chain  models.  These  ideas  are 
then  illustrated  with  simulation  results  and  numerical  calcula- 
tions in  Section  3.  Concluding  remarks  are  made  in  Section  4. 

In  particular  we  discuss  the  possibility  of  combining  internal 
controls  with  external  controls  to  obtain  further  variance  reduction. 


2 . Internal  Controls 
2.1.  General  Ideas 

X = {X^:t  0}  be  the  regenerative  process  being 

simulated.  Recall  that  a process  is  regenerative  if  a renewal 

process  T = ^T^:n  ^ 0}  is  defined  on  the  same  probability  space 

as  X and  if  the  portions  of  the  X process  between  consecutive 

regeneration  points,  the  T 's,  are  i.i.d.  Typically  the  T 's 
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denote  the  times  the  X process  enters  a fixed  state  and  at 
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these  times  the  process  starts  over  from  sci.  ch  independent  <-f 
its  past  history.  For  a formal  definition  of  regenerative  process 
and  general  background  on  the  regenerative  method  of  simulation 
analysis  consult  Crane  and  Iglehart  (1975a) . Let  Tq  = 0 for 
simplicity  and  define  ~ ~ "^k-l'  ^ 2.  ^ portion  of  X 

in  the  interval  t^k-l''^k^  referred  to  as  the  kth  cycle  and 

is  of  length  Tj^  . Suppose  now  that  £{12}  < °°  and  that  the 
distribution  function  of  is  aperiodic  (e.g.  its  support 

is  not  contained  in  a set  of  the  form  {0 ,h, 2h, . . . } , where 
h > 0)  . Then  subject  to  mild  regularity  conditions  =>  X 

as  t Z’  where  =>  means  P{X^  £ x}  ♦-  P{X  £ x}  for  all 

continuity  points  (x)  of  the  limit  distribution.  A similar  result 

holds  when  is  periodic;  see  Crane  and  Iglehart  (19  75a)  . The 

random  variable  X is  frequently  thought  of  as  the  steady-state  con- 
figuration of  the  system  being  simulated.  Suppose  f is  a 

measurable  function  from  the  state  space  of  X to  the  real  line 

and  that  we  are  interested  in  estimating  r = E{f(X)}.  Define 
the  sequence  of  random  variables  (r.v.'s)  {Y,  :k  > 1}  by 


(2.1) 


Y,.  = 


^k 

/ 

'k-1 


f (Xg)ds 


k > 1 . 


If  the  time  parameter  of  X is  discrete  (as  in  a discrete  time 
Markov  chain),  the  integral  in  (2.1)  should  be  replaced  by  a sum. 
The  regenerative  stucture  of  X,  the  process  being  simulated,  gives 
us  the  two  important  properties  stated  in  Section  1 : the  pairs 
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i.i.d.  and  r = E{Yj^}/E{t^}  , provided 

E{|f(X)|}  < oo,  which  we  shall  always  assume.  The  c.l.t.  for  the 

ratio  estimator  r(n)  = Y(n)/T(n)  indicated  in  (1.1)  is  proved 

in  Crane  and  Iglehart  (1975a)  and  follows  from  the  classical 

c.l.t.  for  the  i.i.d.  mean  zero,  finite  variance  r.v.'s 

2 2 

Zj^  = Yj^  - k 1 . We  always  assume  0<a  = E{Zj^}<«>. 

2 

The  variance  of  Z^^,  o , is  also  related  to  the  variance  of  r(n) 
through  the  asymptotic  relation  (as  n ->■  «>) 


(2.2) 


0^{r(n)}  ~ n ^(o^{Z,}/E  {x,})  . 


For  a derivation  of  this  result  see  Cramer  (1964,  p.  354, 
eq.  (27.7.3)).  A number  of  point  estimates  and  confidence  intervals 
have  been  proposed  for  ratios  such  as  r (see  Iglehart,  1975) , but 
the  simplest  conceptually  and  computationally  are  the  so-called 
classical  ones.  The  point  estimate  is  r(n)  and  the  confidence 
interval  is 


(2.3) 


i(n)  = [r(n)  - ^I-y/2^^’^^  ' 


where  z,  = <I>~  (1  - y/2)  , the  inverse  of  the  standard  normal 

1-y/2 

distribution  function,  and  s is  the  classical  point  estimate 
of  a which  is  constructed  as  follows.  Let  s^^^^  [resp.  ^22^ 
be  the  sample  variance  of  the  ® [resp.  Tj^'s]  and  Sj^2 

the  sample  covariance  of  the  ® ^nd  Xj^'s.  Then  s is  defined 

by 

s = [s  - 2(Y/t)s^2  • 


I 
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Introduction  of  the  internal  control  variables  mentioned 

in  Section  1 is  done  with  the  hope  of  reducing  the  variance  of 

the  point  estimator  of  r.  This  in  turn  will  either  reduce  the 

number  of  cycles  that  need  to  be  simulated  for  fixed  precision 

or  reduce  the  length  of  the  confidence  interval,  (2.3),  for  r. 

The  sequence  of  internal  control  variables  1.  defined 

on  the  same  probability  space  supporting  the  process  X in  such 

a way  that  Cj^  only  depends  on  X (or  underlying  r.v.'s 

defining  X)  in  the  kth  cycle.  This  construction  of  the  ' s 

will  insure  that  they  are  i.i.d.  because  of  the  basic  i.i.d. 

structure  of  the  cycles  of  a regenerative  process.  The  control 

variables  allowed,  however,  are  further  restricted  to  those 

for  which  the  mean,  E{Cj^},  is  either  known  or  can  be  calculated 

analytically.  Thus  one  of  the  prices  we  must  be  prepared  to 

pay  to  obtain  variance  reduction  for  our  simulation  is  the 

analytical  work  of  computing  E{C^}  . Having  defined  the 

we  proceed  to  form  new  ratio  estimators  for  r.  Starting  with 

the  ratio  estimator  r(n)  = Y(n)/T(n)  we  can  either 

control  the  Y,  's,  the  t 's,  or  both  if  we  introduce  a second 
^ k 

sequence  of  control  variables.  The  estimator  r^^(n),  proposed 
in  (1.2),  involves  controlling  the  only.  Here  the  sub- 

script CT  is  meant  to  suggest  "con^  ..oiling  the  top."  For 
convenience  let  (for  k 1) 

s'S  - 

and 
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Observe  that  the  Yj^'s  are  i.i.d.  with  E{Yp  = EfY^^}  and  con- 
sequently the  Zj^'s  are  also  i.i.d.  with  E{Zj^}  = 0.  Hence  the 
standard  c.l.t.  yields  (as  n -►  ») 

I zya{Z'}n^/2  = n^/^[  (Y/T)-r]/(a{zM/7)  = N(0,1) 
k=l  ^ -L  1 

which,  upon  replacing  t by  E{tj^}  in  the  denominator  (this 
can  be  justified  by  a continuous  mapping  argument) , oecomes 

n^^^[r  (n)  - r] 

— (a('Z|}A{T^})  ^ N(0,1). 

The  variance  of  Z^,  which  obviously  depends  on  B,  is  easily 
seen  to  be 

a^{Z|}  = o^{Zj^}  + 2B  cov{Zj^,C^}  + B^a^{C^}  . 

2 

Now  select  B so  as  to  minimize  a {Z^}.  This  yields 

cov{Z, ,C, } 

(2.4)  6*  ^ 

a {Cj^} 

and 

(2.5)  a^{Zp  = a^{Z^}  [1  - p^{Z^,C^}]  , 

where  p{z^,C^)  is  the  correlation  coefficient  between  Z^ 

and  C^.  If  we  use  the  c.l.t.  for  the  estimator  r^^(n)  as  a basis  for 

7 


constructing  confidence  intervals  for  r,  then  we  want  a control 

2 

variable  which  will  minimize  a {Z^}.  From  (2.5)  we  see 

2 

that  we  need  a large  value  of  p Since  the  length 

1/2 

of  such  a confidence  interval  is  proportional  to  o{Zj^}/n  ' , 
to  reduce  the  number  of  cycles  required  by  a factor  of  four  we 
need  |p{Zj^,Cj^}|  > (0.75)^^^  = 0.866.  To  obtain  a control,  C^, 
which  is  highly  correlated  with  either  or  is  relatively 

easy  to  do.  However,  because  and  are  themselves  highly 

correlated,  it  is  much  more  difficult  to  find  a control 
which  is  highly  correlated  with  Z^^.  In  general,  we  shall  try 
to  find  a control  which  mimics  Z^  but  for  which  we  can 

still  calculate  its  mean,  eCc^^}.  As  usual  we  try  to  do  as  much 
analytically  as  possible,  before  having  to  resort  to  simulation. 

We  illustrate  some  possible  candidates  for  controls  in  the  context 
of  GI/G/1  queues  and  discrete  time  Markov  chains. 

2.2.  Internal  controls  for  GI/G/1  Queues 

In  the  GI/G/1  queue  we  assume  the  zeroth  customer  arrives 

at  time  t^  = 0 , finds  a free  server,  and  experiences  a service 

time  Vq . The  nth  customer  arrives  at  time  t^  and  experiences 

a service  time  v . Let  the  interarrival  times  t - t . = u , 

n n n-i  n 

n > 1.  Assume  the  two  sequences  (v  :n  > 0}  and  {u  :n  > 1} 

— n — n — 

each  consists  of  independent,  identically  distributed  (i.i.d.) 
random  variables  (r.v.'s)  and  are  themselves  independent.  Let 
E{v^}  = E{u^}  = a and  p = A/p,  where  0 < A,p  < ».  Thus 

8 


W 


y has  the  interpretation  of  the  mean  service  rate  and  X has 

the  interpretation  of  the  mean  interarrival  rate.  The  parameter 

p is  called  the  traffic  intensity  and  is  the  natural  measure 

of  congestion  for  this  system.  We  shall  assume  that  p < 1,  a 

necessary  and  sufficient  condition  for  the  system  to  be  stable. 

While  many  characteristics  of  interest  can  be  estimated 

using  the  regenerative  method,  we  shall  restrict  our  attention 

to  the  waiting  time  of  the  n—  customer,  W (time  from  arrival 

n 

to  commencement  of  service) . For  further  discussion  of  the 
simulation  of  the  GI/G/1  queue  see  Crane  and  Iglehart  (1974b) . 

To  obtain  a representation  for  the  process 

X = V T - u and  set  S,.  = 0,  S =X.  + ***+X,n>l. 
n n-1  n 0 n 1 n — 

The  following  well-known  recursive  relationship  exists  for 

the  W ' s : 
n 


«0  = "„+l  " + =‘n+ll  ' " > “ 


By  induction,  we  also  have 


W^  = max{S^  - Sj^:k  = 0,  1,  ...  , n},  n 0 . 


Using  the  strong  Markov  property  of  the  process  ^S^:n  0} 

it  can  be  shown  that  there  exists  a sequence  of  integer-valued 

r.v.'s  1.  0^  such  that  Tq  = 0,  Tj^  < and  w^  = 0 

k 

with  probability  one.  In  other  words,  the  customers  numbered  T, 


are  those  lucky  fellows  who  arrive  to  find  a free  server  and 
experience  no  waiting  in  the  queue.  The  fact  that  there  exists 
an  infinite  number  of  such  customers  in  the  GI/G/1  queue  is  a 
direct  consequence  of  the  assumption  that  P < 1 and  the 
strong  law  of  large  numbers.  Thus  {W^:n  >^0}  is  a regenerative 
process.  If  we  let  ""  '^k  1'  ^ ^ then  represents 

the  number  of  customers  served  in  the  k—  busy  period  (b.p.) 
and  they  are  numbered  ^"^k-l'  '^k-1  ***  ' define 

the  sequence  {Yj^:k  ^ 1}  by 


\=  I w , 


the  sum  of  the  waiting  times  in  the  k—  b.p.  Since  the  queue 
is  stable  for  p < 1,  we  know  that  =»  W as  n -►  «>  and  wish 

to  estimate  E{W} . 

We  are  interested  in  constructing  controls,  C^,  which  are 
highly  correlated  with  which  are  not  so  complex 

that  we  can  not  calculate  their  means.  The  controls  we  consider 
are  of  the  form  - t^/m,  where  the  r.v.  is  an 

attempt  to  mimic  . To  compute  the  mean  of  we  need,  of 

course,  to  be  able  to  compute  E{Tj^}.  For  M/G/1  queues  we 
know  that  E{tj^}  = l/(l-p).  For  GI/M/1  queues  E{Tj^}  = l/(l-6), 
where  6 is  the  root  inside  the  unit  circle  of  z-E{exp  [-y  (l-z)  Uj^]}  = 0, 
where  u^^  is  an  interarrival  time.  If  the  queue  in  question  is 
neither  an  M/G/1  or  GI/M/1  queue,  then  the  term  T^^/y  in 
will  have  to  be  approximated  by  another  r.v.  whose  mean  can 


be  computed.  The  factor  l/y  multiplying  helps  to  make 

the  variance  reduction  obtained  independent  of  the  scale  parameter 
U.  Also  it  can  be  argued  that  the  term  is  then  in  the 

same  units  as  namely,  the  unit  of  time.  Several  alternatives 

for  are  listed  below,  indexed  by  a superscript.  They  are 


= X-"  - 
- 


Wq  = 0 , 


Wq  + Wi  , 


Ti  = 1 


Ti  >_  2; 


Ti  = 1 


Wq  = 0 , 

Wq  + , 


Wq  + Wi  + X2  , 


1 2; 


^1  = 1 


^1  i 3'- 


°1^^  = + X2)‘^  = ; 


Xi  = 1 


X^  + (xj^  + X^)'^  , 


> 2 


Wq  = 0 , 

Wq  + , 

Wq  + + W2  , 


In  general,  it  is  more  difficult  to  calculate  E{D^^  } as  i 

increases.  On  the  other  hand,  as  i increases  comes 

closer  to  and  presumably  results  in  more  variance  reduction. 

Let  - T^/y,  i = 1,2,3. 

In  our  simulation  runs  for  an  M/M/1  queue  the  controls 

, and  were  generally  found  to  do  much  better  than 

. However,  the  more  complicated  controls,  and 

(2 ) 

gave  little  improvement  over  . Thus  with  both  variance 

(2) 

reduction  and  ease  of  computation  in  mind,  is  our  first 

(2) 

choice.  To  compute  for  the  M/M/1  queue  first  note  that 


Pfx,  = 1}  = P{v^  < u.,  } = / e 
1 u 1 Q 


ye  dy  = (1  + p)  ^ 


E{X^|X^  > 0} 


r°° 

(-^)  J e"^y 

4 + p^  0 


E{xJ^} 


E{Xj|X^  > 0)-  P{Xj  > 0)  = 


Hence 


E{D.^^^}  = O + E{xt  + X^Jlxt  > 0}  = 

1 1+p  1 21 


_P_  [1  + P_ 

1+p  ^y  y(l+p) 


since  X^^  and  are  independent.  The  expectation  of  C^‘ 

is  then  given  by 
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Further  discussion  on  the  computation  of  the  mean  of  controls  such 
as  for  the  GI/G/1  queue  is  contained  in  the  Appendix.  Com- 

putational results  from  our  simulation  runs  are  contained  in 
Section  3. 


2.3.  Internal  controls  for  Markov  chains 

The  second  example  we  consider  is  a discrete  time  Markov 
chain.  Assume  now  that  we  are  simulating  an  irreducible,  aperiodic, 
positive  recurrent  Markov  chain,  X = {x^:n  ^ 0},  with  the  goal 
of  estimating  r = E{f(X)}.  The  controls  we  consider  here  are 
of  the  form 


noA(Ti-l) 


^1  = 


[f (Xn ) - r.]  , 


where  r^  is  a guess  of  r,  n^  a fixed  integer,  and 

Hq  a (t^-I)  denotes  the  minimum  of  n^  and  (t^-1) . Again  our 

motivation  for  is  to  mimic 

2 = I t£(X  ) - rl  . 

«,=0 

Of  course,  we  still  must  be  able  to  compute  E{C^}  in  order 
to  implement  the  method  of  internal  control  variables.  Let  ^ 
be  the  state  space  of  X,  which  we  assume  is  finite,  and  P the 


transition  probability  matrix.  Furthermore,  let  g be  the 
column  vector  with  components  f(i)  - r^  and  be  the  matrix 

with  elements 


I Pk«-' 
jPkJ,  I 

' 0 , 


Z ^ j 
£ = j . 


Then  if  we  let  = j and  form  regenerative  cycles  based  on 
the  times  of  return  to  state  j,  it  can  easily  be  shown, using  the 
method  developed  in  Hordijk,  Iglehart,  and  Schassberger  (1976), 
that 


(2.7) 


E{C,}  = { I 
^ n=0 


=n  J J 


where  j is  the  regenerative  state  and  the  subscript  j means  the 
jth  component  of  the  indicated  vector.  For  two  simple  Markov  chains 
the  theoretical  amount  of  variance  reduction  obtained  using  the 
control  has  been  calculated  for  different  regenerative  states  j 

and  values  n^  and  is  contained  in  the  next  section.  From  (2.6) 
it  is  clear  that  the  larger  the  integer  n^  selected  the  closer 
comes  to  duplicating  7,^.  However,  the  larger  nQ  the  more 
difficult  is  the  computation  of  E{C^}  contained  in  (2.7). 
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3.  Numerical  Results 


3.1.  The  M/M/1  queue 

Even  for  the  simple  M/M/1  queue  it  is  in  general  impossible 
to  verify  analytically  what  variance  reduction  will  be  obtained 
via  the  several  internal  controls  listed  in  the  previous  section, 
or  to  get  an  idea  of  the  magnitude  of  the  effect.  For  something 
as  simple  as  it  is  difficult  to  compute  analytically  the 

correlation  between  and  for  the  M/M/1  queue,  and  this 

is  what  is  required  in  equation  (2.5)  to  find  the  variance 
reduction. 

Thus,  we  resorted  to  simulations  to  verify  the  amount  of 
variance  reduction  obtained  and  the  relative  effectiveness  of 
the  various  controls.  In  the  final  simulations  all  runs  were 
performed  on  an  IBM  System  360/67  computer  using  the  LLRANDOM 
package  (Learmonth  and  Lewis  (1973))  which  generates  random  numbers 
according  to  the  scheme  given  by  Lewis,  Goodman,  and  Miller  (1969) 
and  exponentially  distributed  random  numbers  using  the  Marsaglia 
"rectangle-wedge- tail"  method.  Tests  of  the  random  number 
generator  are  given  in  Learmonth  and  Lewis  (1974) . 

Of  the  extensive  simulation  runs  performed,  we  give  here 
only  a summary  of  the  conclusions  and  two  detailed  tabulations 
in  the  case  of  the  most  suitable  control. 

(1)  The  controls  » ^1^^  ^1^^  much  better  generally 

than  with  little  improvement  over  obtained  by 

use  of  and  We  say  generally  because  results 

vary  unpredictably  with  X and  y and  their  ratio  p. 
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(2)  Subtracting  the  number  of  customers  served  in  a busy  period 

generally  improves  the  variance  reduction.  By  making  it 

dimensionally  stable  as  in  we  obtain  a "variance 

reduction"  measured  in  terms  of  ratios  of  standard  deviations, 

of  approximately  70%,  uniformly  over  X and  y.  This  is 

roughly  equivalent  to  halving  the  number  of  cycles  (b.p.'s) 

2 

that  one  must  simulate;  (0.7)  ~ .5.  Much  better  reductions 

can  be  obtained  for  smaller  p (i.e.  p = 0.25)  by  specially 
designed  controls;  the  point  is  that  works  even  out 

at  p = 0.99,  where  variance  reduction  is  extremely  important. 

Table  1 shows  results  obtained  by  simulating  an  M/M/1 

(2) 

queue  out  to  n = 2000  cycles  with  control  and  replicating 

the  simulation  either  m = 250  or  m = 100  times  to  estimate  the 
variance  of  the  estimators  r(n),  r^^(n),  and  r^g(n),  where  we 
drop  the  n for  convenience.  Here,  we  have  specifically  that 


rcB(n)  = 


1 « 

n J/lc 


E(c/’  ))  1 


The  estimated  precision  (standard  deviations)  of  the  estimates 
of  E(W)  are  given  in  brackets  under  the  estimates. 

The  results  in  Table  1 are  for  p = 0.5  (and  m = 250) 
and  two  values  of  u and  for  p = 0.99  (and  m = 100) 
for  two  values  of  y commensurate  with  those  given  for 
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p = 0.5.  The  second,  third  and  fourth  columns  in  the  Table  give 

estimates  of  the  correlations  between  the  control  and  etc.  from 

which  the  theoretical  variance  reduction  can  be  computed.  They 

are  very  close  to  the  values  given  in  the  next  to  last  column, 

the  observed  variance  reduction,  from  which  we  deduce  that 

estimating  3,p  and  Bg  affects  the  variance  reduction  only 

slightly.  Overall  there  is  negligible  effect  of  different  values 

of  y on  the  variance  reduction  obtained  for  the  cases  p = 0.5 

and  p = 0.99.  For  the  results  p = 0.99  given  in  Table  1,  the 

variance  reduction  is  75%,  which  is  about  the  same  as  for 

p = 0.5.  For  the  case  where  the  control  is  on  the  bottom,  i.e. 

for  T,  the  variance  reduction  is  not  quite  as  good  for  control 

of  Y.  Note  too  that  the  estimated  values  of  E{W}  appear  in 

some  cases  to  be  at  least  three  or  four  standard  deviations  from 

the  true  mean.  This  is  because  the  estimates  r,  r and  r 

Cl  Co 

can  be  seen  from  the  100  replications  to  be  non-normal.  In  other 
words,  for  high  p(0.99),  the  simulation  needs  to  be  taken  out 
further  than  2000  cycles.  In  summary,  the  variance  reduction 
obtained  with  the  internal  control  is  practically  independent 
of  the  scale  factor  y and  the  traffic  intensity  p. 

3.2.  Two  simple  Markov  chains 

We  turn  now  to  two  simple  Markov  chains  to  give  another 
illustration  of  the  method  of  internal  controls.  The  first 
chain  is  simply  a recurrent  random  walk  with  state  space 
E = {0,1, 2, 3, 4}  and  transition  matrix 
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Suppose  we  are  simulating  to  estimate  r = E{ X} , which  can  be 

computed  to  have  the  value  2.  Table  2 contains  the  theoretical 
2 2 

values  of  a {Z^}/a  {Z^}  and  p(C^,Z^),  where  is  given  by 

(2.6)  for  several  choices  of  np  and  rQ . The  entries  in  these 
tables  were  calculated  using  the  methods  in  Hordijk,  Iglehart, 
and  Schassberger  (1976) . Note  that  the  variance  reduction  is 
largest  when  E{t2^}  is  smallest  and  np  is  largest.  This  also 
yields  the  highest  value  of  p(C^,Zj^)  in  )ceeping  with  our  idea 
of  having  mimic  . Not  too  much  is  lost  in  variance 

reduction  as  the  guess,  r^,  departs  from  the  true  value  r = 2. 
A variance  reduction  of  more  than  50%  is  attainable,  but  depends 
on  the  state  which  is  used  as  the  regeneration  point. 

The  second  Markov  chain  represents  an  (s,S)  inventory 
model;  see  Crane  and  Iglehart  (1974)  for  further  details.  The 


= {6,7, 

. . . ,10} 

and  the 

P 

matrix 

3/8 

0 

0 

0 

5/8 

1/4 

3/8 

0 

0 

3/8 

3/16 

1/4 

3/8 

0 

3/16 

1/8 

3/16 

1/4 

3/8 

1/16 

1/16 

1/8 

3/16 

1/4 

3/8 

71.474  49831.35  81%  1000.00 

22.323) 


r 


TABLE  2 
2 2 

Variance  ratio  a /a  {Z.^}  (top  row) 


and 

Correlation 

p(C^,Zl) 

(bottom  row) 

for  the 

random  walk 

model 

$ 

ro  = 1.5 

Return  State 

i E{t^} 

0 

8 

0.97 

1.0 

0 .91 

-0.16 

-0.01 

0.30 

1 

4 

0.82 

0.73 

0.57 

0.43 

0.52 

0.65 

2 

4 

0.62 

0.51 

0.40 

0.62 

0.70 

0.78 

3 

4 

0.94 

0.91 

0.84 

0.25 

0.31 

0.40 

4 

8 

0.97 

0.98 

0.99 

-0.16 

-0.16 

-0.10 

"0  = 2 

0 

8 

0.97 

0.98 

1.0 

-0.16 

-0.12 

0.03 

1 

4 

0.89 

0.83 

0,70 

0.34 

0.42 

0.55 

2 

4 

0.57 

0.46 

0.35 

0.66 

0.74 

0.81 

3 

4 

0.89 

0 .83 

0.70 

0.34 

0.42 

0.55 

4 

8 

0.97 

0.98 

1.0 

-0,16 

-0.12 

0.03 

TABLE  3 


2 2 

Variance  ratio  a {Z^}/a  { Z^}  (top  row) 
and  Correlation  p(Cj^,Z^)  (bottom  row) 
for  (s,S)  Inventory  Model 


o 

II 

00 

Return  State  i 

^ 2 

"0=2 

no=4 

6 

5.55 

0.75 

0.65 

0 .49 

0.51 

0.59 

0.72 

7 

5.72 

0.82 

0 .67 

0.49 

0.43 

0.57 

0.72 

8 

6.27 

0.89 

0.76 

0.57 

0.33 

0.49 

0.66 

9 

7.21 

0.95 

0.86 

0.69 

0.21 

0.38 

0 .56 

10 

2.88 

0.95 

0.81 

0.42 

0.22 

0.44 

0.77 

o 

II 

6 

5.55 

0.75 

0.75 

0.78 

0.51 

0.50 

0.47 

7 

5.72 

0.88 

0 . 80 

0.73 

0.34 

0.45 

0.52 

8 

6.27 

0.93 

0.83 

0.71 

0.27 

0.41 

0.54 

9 

7.21 

0.95 

0.86 

0.72 

0.23 

0.37 

0.53 

10 

2.88 

0.79 

0.57 

0.24 

0.46 

0.66 

0.87 
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Again  we  wish  to  estimate  E{X}  = 8.297.  Table  3 contains  the 
results  analogous  to  those  for  the  random  walk  (Table  2) . The 
same  general  observations  about  and  n^  made  for 

Table  2 seem  to  hold  here  as  well.  Note  however  that  the 
results  are  sensitive  to  the  choice  of  Tq . Again  a 50% 
reduction  in  variance  can  be  attained. 

4 . Conclusions  and  Extensions 

We  have  been  able  to  obtain  a 50%  variance  reduction  using 
internal  control  variables,  for  the  regenerative  estimate  of  the 
limiting  value  of  the  mean  waiting  time  in  an  M/M/1  queue.  This 
reduction  is  obtained  uniformly  over  all  parameter  values.  It  is 
fairly  certain  that  the  technique  will  work  well  with  any  GI/G/1 
queue . 

Internal  control  variables  can  be  easily  used  with 
discrete  time  Markov  chains.  The  examples  used  in  this  paper 
showed  that  a variance  reduction  of  50%  is  attainable.  This  figure 
is  likely  to  vary  widely  with  the  particular  Markov  chain. 
Continuous  time  Markov  chains  and  semi-Markov  processes  can  be 
handled  in  the  same  way  using  the  discrete  time  method  of 
Lordijk,  Iglehart  and  Schassberger  (1976)  . 

Another  method  of  internal  stratified  sampling  was  also 
investigated.  This  method  produced  little  overall  variance 
reduction  despite  considerable  effort. 
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Finally  we  note  that  it  is  possible  to  apply  the  idea  of 
internal  controls  to  the  classical  sample  average  estimate  over 
a realization  of  fixed  length  m.  Thus  in  estimating  E(W) 
in  the  GI/G/1  queue  we  have 


, m-1 

W(m)  - - y W.  , 

TTI  - “ _ *1 


(4.1) 


which  may  be  written  as 


I N(m) 

W(m)  = i { I Y.  + Y 


n(m)+l^  ' 


(4.2) 


where  N (m)  is  the  number  of  completed  busy  periods  in  the  queue 
in  [0,m],  Yj^  as  before  is  the  sum  of  the  waiting  times  in  the 
Icth  cycle,  and  ^N(m)+1  waiting  times  in  the 

last,  incomplete  cycle. 

A central  limit  theorem  for  W{m)  holds  for  which  the  variance 

2 

term  is  proportional  to  E{Zj^},  and  which  is  now  estimated 
from  the  random  number  N (m)  of  cycles.  Control  is  applied 
to  the  Yj^'s  in  (4.2)  just  as  it  is  applied  in  the  ratio 
estimator  r^^(n).  Call  this  estimator  (m) . For  the  M/M/1 
queue  the  variance  reduction  observed  in  the  simulations  was 


the  same  for  all  controls  C^ 
with  (m) . 


with  the  ratio  estimator  and 


The  main  reason  for  considering  W^(m)  is  that,  while 
it  loses  the  advantage  of  being  an  estimator  using  a fixed  number 
of  i.i.d.  random  variables,  one  can  apply  the  classical  external 
controls  to  W^(m).  Thus  one  could  use  the  difference  of  the  sum 
of  the  m arrival  times  and  the  m service  times  to  control 
W^(m).  We  have  not  yet  tried  this  idea. 
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To  implement  the  internal  control  techniques  for  a 
GI/G/1  queue,  certain  theoretical  parameters  are  required.  In 
this  appendix  we  shall  indicate  the  values  of  these  parameters 
in  so  far  as  they  can  be  calculated.  These  values  are  either 
well-known  or  easily  calculated.  For  a reference  to  the  known 
formulas  see  Cohen  (1969)  . 

We  begin  with  E{t^},  the  expected  number  of  customers 

served  in  a busy  period.  For  the  general  GI/G/1  queue  recall 

that  we  let  X = v n - u and  S = X,  + • • • + X , for  n > 1, 
n n-1  n n 1 n — ' 

with  S„  = 0 . Then  t,  = inf{n  > 0:S  < 0}.  The  general 

0 1 n — 

expression  for  E{tj^}  is  given  by 

oo 

E{t,}  = exp{  y n~^P[S  > 0]}  , 

1 ^ ^ n 

n=l 

an  impossible  expression  to  evaluate  in  general.  Another  useful 
expression  for  E{tj^}  is 

(A.l)  E{t^}  = 1/P{W  = 0}  , 

where  W is  the  stationary  waiting  time.  In  the  special  case 
of  the  M/G/1  queue,  however,  we  have 

E{t^}  = (1  - p)'^  . 


•J. 
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Now  for  the  queue  GI/M/1  we  can  use  (A.l)  and  the  stationary 
distribution  of  the  embedded  Markov  chain  to  conclude  that 


= (1  - 6) 


-1 


where  6 is  the  root  inside  the  unit  circle  of 

z - U{u (1-z) } = 0 


-su 


with  U(s)  = E{e  Re  s ^ 0,  and  where  Vg  is  an  exponential (y) 

r,v.  It  is  easy  to  check  that  6 = p for  M/M/1  queues.  Daley 
(1975)  has  recently  proposed  the  approximation  to  6 given  by 


-1, 


-1 


6 = a^(l-p)"  + 2(l-b  ■")p  + (2b  -1)  p , 

where  a^  = P{Uj^=0},  E{u^}  = 1,  and  b = E{u^}  . This  approximation 
gives  good  results  in  a number  of  examples  calculated  by  Daley 
(1975)  and  may  be  useful  for  the  purposes  of  this  paper. 

Next  we  turn  to  the  computation  of  P{tj^=1}  and  P{tj^=2} 
For  the  GI/G/1  case  we  have 


P{tj^=1}  = P{S^  < 0} 


and 


P{t^=2)  = P{S^  > S2  1 0}, 


25 


both  of  which  can  be  worked  out  with  a little  effort.  For  the 
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M/M/1  queue 

P{t^=1}  = (1+p)'^ 

and 

P{tj^=2}  = p(l+p)  ^ . 

FoJ  the  M/G/1  queue 

P{Tj^=1}  = V(A)  , 

-XV 

where  V(X)  = E(e  ),  while  for  the  GI/M/I  queue 

P{t^=1}  = 1 = U(y)  , 

where  U(s)  is  given  above.  For  the  M/Ej^/1  and  Ej^/M/1  queues 

the  value  P{Tj^=2}  can  be  calculated  with  some  effort.  As  these 

expressions  are  cumbersome  they  shall  be  omitted. 

Finally  we  give  various  partial  expectations  which  are 

(4 ) 

needed  for  computing  the  means  of  internal  control  . Namely, 

E{S^  + S2,  I 2}  = E{Sj^,  > 0}  + E{S2,  > 0} 

and 

E{Ti,  >_  2}  = E{Tj^}  - P{Tj^=1}. 

Here  the  symbol  E{X,A}  = E{Xl^} , where  X is  a r.v.,  A an 
event,  and  1.  the  indicator  function  of  A.  In  the  special 
case  of  the  M/M/1  queue. 
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